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Systems of charged particles on anisotropic three-dimensional lattices are investi- 
gated theoretically using Debye-Hiickel theory. It is found that the thermodynamics 
of these systems strongly depends on the degree of anisotropy. For weakly anisotropic 
simple cubic lattices, the results indicate the existence of order-disorder phase transi- 
tions and a tricritical point, while the possibility of low-density gas-liquid coexistence 
is suppressed. For strongly anisotropic lattices this picture changes dramatically: the 
low-density gas-liquid phase separation reappears and the phase diagram exhibits 
critical, tricritical and triple points. For body-centered lattices, the low-density gas- 
liquid phase coexistence is suppressed for all degrees of anisotropy. These results 
show that the effect of anisotropy in lattice models of electrolytes amounts to reduc- 
tion of spatial dimensionality. 
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I. INTRODUCTION 

Understanding thermodynamic properties of electrolyte systems is a long-standing 
problem!!'! which, in recent years, has attracted increased attention due to controversial 
results on the nature of criticality in Coulomb systems.!'! ExperimentsB suggest that the 
critical region of solutions of some organic salts can be well described by classical behavior, 
while for other electrolyte systems the Ising-like description with non-classical behavior is 
more appropriate.! This dichotomy has greatly stimulated theoretical attempts to under- 
stand the thermodynamics of ionic systems. 

Criticality in simple nonionic fluids can be successfully described and analyzed by the 
renormalization group (RG) method, and it is reasonable to suggest that this method could 
also be used to investigate also Coulomb systems. However, in order to proceed with RG 
calculations, a physically meaningful and well based mean-field theory should be developed.! 
Most theoretical studies of ionic systems concentrate on the simplest model of electrolytes, 
the so called restricted primitive model (RPM), in which ions are viewed as equal size 
particles of positive and negative charges of equal magnitude. Currently, there are two the- 
oretical directions in the development of the mean-field description of ionic fluids. The first 
approach is based on integral equations for correlation functions!ffl!'0, while the second ap- 
proach extends the original Debye-Huckel (DH) theory.i'SIillillil Comprehensive theoretical 
analysis, which utilizes for example thermodynamic energy bounds, and comparison 
with current Monte Carlo simulations,0'0 indicate that theories based on DH theory may 
provide a better description of the thermodynamics of electrolytes in critical regions. 

So far most of the theoretical efforts in the investigation of charged systems have been 
devoted to continuum models. However, lattice models are also important for understanding 
the criticality in Coulomb systems, since the Ising model, which is a lattice gas model, 
has been crucial for the description of critical phenomena in nonionic systems.!'! There 
are few numerical0'0 and analyticallBQ results for lattice ionic systems which show that 
the phase diagram differs significantly from continuum models. The structures of simple 
cubic (sc) and body-centered cubic (bcc) lattices allow for charge distribution with the 
appearance of a long-range order phase at low temperatures similar to that of an ionic 
crystal. However, this ordering decreases entropy and for high temperatures the disordered 
phase is thermodynamically more stable. As a result, there is an order-disorder phase 
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transition line which ends up at the tricritical point, EHifl while for continuum RPM 
systems only a gas-liquid coexistence can be found. 

Recently, a systematic investigation of lattice models of electrolytes has been presented.0 
In this investigation the lattice restricted primitive model (LRPM), with charged particles 
occupying sites on a general d- dimensional lattice, has been considered using Debye-Hiickel 
theory. By solving exactly the lattice version of the Debye-Hiickel equation, closed ex- 
pressions for thermodynamic properties of general d- dimensional ionic systems have been 
obtained. For three-dimensional lattice Coulombic systems specific calculations, which in- 
cluded pairing and dipole-ion solvation, yielded a gas-liquid phase separation at low densities. 
However, by taking into account the lattice symmetry, it has been shown that for sc and bcc 
lattices this gas-liquid phase separation is thermodynamically unfavorable, and the order- 
disorder phase transitions with the tricritical point will dominate, in agreement with Monte 
Carlo simulation results.0lll 

The thermodynamics of lattice anisotropic ionic systems have not been studied yet, al- 
though they may provide important information on thermodynamics of real electrolytes. 
In addition, the lattice stretching, which leads to anisotropy, can be viewed as analogous 
of lowering the spatial dimensionality of the system. However, DH-based calculations for 
the continuumS and latticed models of electrolytes predict an increase of gas-liquid critical 
temperatures for lower dimensions. Thus, for lattice anisotropic models, gas-liquid phase 
separation may reappear along with distinct order-disorder phase transitions. This possibil- 
ity raises the question of the precise determination of phase diagrams for ionic systems on 
anisotropic lattices. 

In this article, we present a theoretical investigation of anisotropic lattice models of 
electrolytes using the Debye-Hiickel method,@ which has the advantage that it accounts for 
both electrostatic screening and sublattice ordering in a unified framework. We investigate 
the ionic systems on three-dimensional lattices obtained by anisotropic stretching of simple 
cubic and body-centered cubic lattices. The paper is organized as follows. In Sec. II we 
provide a combined DH and mean-field ordering description of ions on simple tetragonal 
lattices. The similar analysis for stretched body-centered lattices is given in Sec. III. A 
discussion and conclusions are presented in Sees. IV and V. 
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II. DEBYE-HUCKEL THEORY FOR SIMPLE TETRAGONAL LATTICES 

A. Pure DH theory 

Our derivation for anisotropic lattice electrolytes follows closely the Debye-Hiickel ap- 
proach for isotropic lattice ionic systems.0 We consider a system of equal numbers of pos- 
itive and negative ions with the total density p = p_ + p + on a tetragonal lattice, which is 
obtained by stretching the simple cubic lattice, with unit cell dimensions a x a x b with an 
anisotropy factor defined as the ratio of lattice parameters, a = a/b. The case a = 1 cor- 
responds to the isotropic simple cubic lattice electrolytes which have been studied in detail 
earlier.0 The Debye-Hiickel approach implies that we can construct the total free energy 
of the system by summing consecutively the terms which describe interactions between dif- 
ferent species. Charged particles interact through the lattice Coulomb potential, otherwise 
they behave as ideal particles with additional hard-core on-site exclusions. Thus the total 
free energy density is given by / = f Id + f DH - The ideal lattice gas contribution can be 
written as 

/" = -7^7 = -- " Ml - Pi, (1) 

k B TV v v 

where p* = pv is the reduced dimensionless density and v = a 2 b is the unit lattice cell 
volume. 

The other term in / comes from the Coulombic interactions of the free ions and includes 
the effect of screening. In order to find an expression for the Debye-Hiickel contribution to 
the free energy density, we need the electric potential felt by an ion due to all other ions. 
This potential can be found by fixing an arbitrary ion at the origin and solving the linearized 
Poisson-Boltzmann equation 

A^(r) = «V(r) - (C d q/Dv )5(r), (2) 

where k 2 = Cdl3pq 2 /D is the inverse squared Debye screening length, with (3 = 

and Cd = 47r for three-dimensional lattices. The lattice Laplacian used in Eq.@ can be 

presented in the form which incorporates the geometry of the lattice, 

Ay? = A x ip + A y Lp + A z ip, (3) 
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with 

A^(r) = 1/a 2 [<p(r - a^) - 2<p(r) + cp(r + a^)} , (4) 

where % 6 and are the unit vectors along the corresponding 

lattice directions. Then Eq.(0) can be easily solved by Fourier transformation to yield for 
r > 

^ = 3Z^ J (x 2 + 4 + 2« 2 )/6- J(k,a 2 )' (5) 
where we introduced the anisotropic lattice function 

J(k, 7) = — V e ik a = - (cos fa + cos fc 2 + 7 cos fc 3 ) , (6) 



mi 



Co = 6 is the number of nearest neighbors, and j k = (27i)~ d d k and 2 = «a. Note that 
the vectors k = {fa/a, kija, k^/b) describe the reciprocal lattice. 

To compute the electric potential at the origin due to the surrounding ions, it is recalled 
that no other ion can be placed at that site, and hence the appropriate equation is the 
Laplace's equation, 

A^(r = 0) = 0. (7) 

This enables us to write 

A^(r = 0) = A^(r = 0) = A^(r = 0) = 0, (8) 



which yields 



ip(0) = (ip(r - ae x ) + <^(r + ae x ))/2 

= (ip{r-ae y ) + tp(r + ae y ))/2 (9) 
= ((p(r - ae z ) + (p(v + ae z ))/2 



Then from Eqs. (El) and (0) it follows that 



2nqa 2 f J(k, a 2 ) 

^ ' = Dv (2 + a 2 ) J (x 2 + A + 2a 2 )/Q- J(k, a 2 ) ' ' ° ! 



k 
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Following the DH approach for isotropic lattices,^ we introduce the integrated anisotropic 
lattice Green's function 

1 



P(z,a) 



1 — z J(k, a 2 



(11) 



which has been evaluated exactly in terms of a product of two complete elliptic integrals of 
the first kind by Joyce.El Then the total potential at the origin due to the surrounding ions, 
ip(0) = ip(0) — f(0)\ x=0 , takes the form 



1> 



2nq 



Db(2 + a 2 



P 
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a) -P 
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a 



(12) 



x 2 + 4 + 2a 2 ' J \4 + 2a 2 ' 

By using the Debye charging procedure, the reduced electrostatic free energy density can 
be calculated explicitly, yielding 



DH 



1 



4(2 + a 2 )v 



x 2 P 
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4 + a 2 



a 



P 
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x 2 + 4 + a 2 



, a I d(x 2 ) 



(13) 



The chemical potential, /x = [i/k B T = —df /dp, is then given by 



/z = lnp*-ln(l-p*)- 



7T 



P 



{2 + a 2 )T* [ \A + a 2 ' 



6 



a)-P 



x 2 + 4 + a' 



(14) 



where, following the continuum DH theory@ and the DH theory for isotropic cubic lattices,0 
the reduced temperature is defined as 

Dk B Tv Db Da 



q 2 a 2 q 2 /3 q 2 (3a ' 



and for the reduced density we obtain 



x 2 T* 
Air 



(15) 



(16) 



Knowing the free energy density and the chemical potential allows us to calculate the pres- 
sure, p = p/ksT = max p [/-|- pp], yielding 



pv 



ln(l-p*) + 



4(2 + a 2 



x 2 P 



6 



x 2 + 4 + 2a- 



■, a 



P 



6 



x 2 + 4 + 2a 2 



, a \ d{x 2 ) 



(17) 



Eqs. (P, (|i~3|), (|14D and ( p!7[) provide a full thermodynamic description of the simple 
tetragonal lattice model of electrolytes. The thermodynamics at the critical region can be 
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investigated by analyzing the spinodal, which is determined by the condition p^ = 0. Using 
Eq. fll^D we obtain 

s (2 + a 2 ) 2 + (l-C) 2 9P(C)/aC 1 ] 

with C = 6/(x 2 + 4 + 2a 2 ). 

The phase transitions and the gas-liquid coexistence can be studied by analyzing the 
pressure and the chemical potential in different phases. The predicted gas-liquid coexistence 
curves for simple tetragonal lattices are shown in Fig.l. The critical temperature increases 
monotonically as the anisotropy parameter decreases and reaches the value of T* = 1/2 at 
a = 0: see Fig. 2a. At the same time, the critical density shows a non-monotonic behavior 
with two maxima and a minimum, and finally approaches the value p* = at a = 0, as 
shown in Fig.2b. Lowering the anisotropy parameter can be visualized as the stretching 
the lattice along one direction, and the limit of a — > corresponds to an infinite distance 
between the layers. Thus the anisotropic lattice model at a = is equivalent to the 2D 
Coulomb system on square lattice, for which the pure DH theoryQ predicts the critical 
parameters to be T* = 1/4 and p* = 0. The apparent discrepancy between our results 
for the critical temperature and the results for two-dimensional lattice electrolytes^ can be 
easily explained by analyzing the linearized Poisson-Boltzmann Eq.(^). In our calculations 
we used the three-dimensional coefficient Cd = 47r, while in two dimensions this coefficient 
is equal to 2tt, which explains the factor 2 in the difference in the corresponding values of 
the critical temperatures. 

At the limit of large a the lattice is stretched along 2 directions, and thus the anisotropic 
lattice with a = oo corresponds to the ID lattice Coulomb system. Our predictions for 
critical parameters in this case are T* = and p* = 0, while for one-dimensional lattice 
electrolytes the DH-based calculations give T* = oo and p* = 0. In this case, the difference 
in critical temperatures can be attributed to our definition of the reduced temperature in 
Eq. (|T5D . Note, however, that the DH method is incorrect in describing one-dimensional ionic 
systems.0 The overall agreement between our estimates of critical parameters of strongly 
anisotropic lattice models of electrolytes and the results for ID and 2D ionic lattice sys- 
tems supports our arguments that anisotropic stretching is analogous to lowering of spatial 
dimensionality for simple cubic lattices. 

A surprising feature of the predicted coexistence curves is the critical density dependence 
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on the lattice anisotropy as exhibited in Fig. 2b. It shows two maxima (at a = 0.4 and 
a = 4.175) and one minimum (at a = 1); this picture is probably the result of geometric 
and packing effects. 

B. Sublattice Ordering 

The Debye-Hiickel approach presented above describes only the low-density behavior of 
the system. To obtain the full thermodynamic description of simple tetragonal lattices we 
have to take into account the lattice symmetry. A simple tetragonal lattice, similarly to 
a simple cubic lattice, can be viewed as consisting of two intercalated sublattices. Ions of 
opposite signs can be distributed unequally between these sublattices, thus reducing the 
electrostatic contribution to the free energy. At the same time, unequal distribution of 
charged particles lowers the entropy which leads to an increase in the total free energy. The 
competition between these factors determine the thermodynamics and phase behavior of the 
system. 

We consider again a simple tetragonal cubic lattice with N charged particles. The overall 
system is neutral, and there are Nj[ (iVj) positive (negative) particles in sublattice A, and 
Ng (Ng) positive (negative) particles in sublattice B. Assuming that sublattice A has an 
excess of positive ions, the corresponding order parameter can be defined as 

N+ + N A N+ + N B - 

This order parameter has a positive value in the ordered phase, while it equals to zero in 
the disordered phase. 

The nonzero charge density on each sublattice produces an additional "background" 
potential $(r) which, however, does not change the correlation functions.0 This potential 
can be found using the linearized Poisson-Boltzmann equation 

A$(r A ) = -{4ir/D)pyq. (20) 

Because of the symmetry between sublattices, we have <&(ta) = — $(rg). By using the 
definition of the lattice Laplacian ©-(§]), and following the approach outlined for isotropic 
lattice electrolytes,0 we obtain for the potential tp due to all ions except the one fixed at 
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the origin the following expression 

- -ihrm + ^" < 21 > 

where ip DH is given in Eq. (|T2|) . The electrostatic part of the total free energy then follows 
again from the Debye charging process, while the entropic contribution can be calculated as 
for isotropic lattices,0lil yielding / = f Id + f DH + f° rd with 

f° rd = * 2 /X* - f- K 1 + ^ ^ + y) + 0--y) ln d - y) - 2in2] , (22) 

2(2 + cr) v T* 2vq 
where f Id and f DH are given by Eqs.flD) and fll3j). 

The knowledge of the total free energy allows us to investigate the possibility of sublattice 
ordering. It can be done by looking for minima of f° rd for nonzero values of the order 
parameter y. This procedure leads us to the equation describing the A-line, along which 
second-order phase transitions occur, 

2 + a 2 

p* = — — — T*. (23) 

7T 

The anticipated tricritical point can be found by calculating the intersection of the A-line 
with the spinodal dp/ dp* = 0, and this analysis yields the equation for the tricritical 
point ,0-0 



4(2 + a s 

Ptri 



8P s [6/(x 2 + 4 + 2a 2 ), a] 



+ = 0. (24) 

1 — n* 9 

2 =4(2+a 2 ) 1 Ptri Z 



d(x 2 ) 

The resulting tricritical densities and temperatures for different values of the anisotropic 
parameter a are presented in Figs. 2a and 2c. The tricritical temperature is a decreasing 
function of the anisotropy parameter, and it reaches its maximal value of T tri = 0.5960 at 
a = 0, which is exactly twice the value of the tricritical temperature for the ionic system on 
the two-dimensional square lattice (see Eq.(70) of Ref.0). This deviation is again the result 
of using different dimension-dependent coefficients Cd, as was argued above. Both critical 
and tricritical temperatures vanish at large anisotropies, however, T tr i becomes smaller than 
T c for a > 4.25. The behavior of the tricritical density is different. It has a minimal value 
for the isotropic lattice (a = 1), and it reaches the maximal values of p c = 0.3794 and 
p c = 0.416 for a = and a = oo respectively. 

The phase diagrams for simple tetragonal lattices are presented in Fig. 3. For weakly 
anisotropic lattices there are only order-disorder phase transitions, while for strongly 
anisotropic lattices the gas-liquid phase separation reappears at low densities. 



10 



III. THEORY FOR BODY-CENTERED TETRAGONAL LATTICE 



Consider a system of equal numbers of positive and negative ions on the body-centered 
tetragonal lattice with 2a x 2a x 2b unit lattice cell. Using the symmetry of the lattice and 
applying Eqs.(|^) and the lattice Laplacian for tetragonal body-centered lattice is given 
by 

2(2 + a 2 ) ■ 



A<p(r) 



3c a 2 



[(p(r + & nn ) - ip(r)\ , 



(25) 



where a = a/b and the summation runs over all Cq = 8 neighbors in a cell. When a = b this 
reduces to the well known lattice Laplacian for the body centered cubic lattice.0@ Because 
of the special symmetry of body-centered lattice, the lattice function 

J fc (k) = i-X\ 



ik-a 



cos ki cos k 2 cos k 3 , 



(26) 



in contrast to the simple tetragonal lattice, is independent of the anisotropy parameter a. 
Thus the isotropic bcc lattice Green's function Pb{z) = f k 1 _ z j b ^ can be used for calculation 
of the thermodynamic properties. Then the potential at the origin due to the surrounding 
ions takes the form 



2nq 3a 2 
Dv (2 + a 2 



Ph 



4 + 2a 2 



x 2 + 4 + 2a 2 



and the electrostatic free energy is given by 

1 



/ 



DH 



4(2 + a 2 )v 



x 2 P h (l) 



P, 



4 + 2a 2 



d{x 2 



(27) 



(28) 



x 2 + 4 + 2a 2 , 

Furthermore, the ordering free energy does not depend on the type of lattice0 and is the 
same both for the simple and body centered lattices. Therefore we can use the corresponding 
expression from (^2|) and estimate the trier itical point from Eq.(p4]) by using Pb{ x i+l+l a s ) 
instead of the simple tetragonal Green's function P[6/(x 2 + 4 + 2a 2 ), a]. Phase diagrams 
for electrolytes on the tetragonal body-centered lattices are shown in Fig. 4. 

Analysis of phase diagrams for the body-centered lattice models of ionic systems indicate 
that both critical and tricritical densities are independent of the degree of lattice stretching, 
while T* and T^ ri are decreasing functions of a. These relations can be understood by 
analyzing the corresponding equation for the spinodal, 

2tt C(l - OdP b (()/d( 



(2 + a 2 ) 2 + (l-C) 2 dP b (0/dC 



(29) 
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where ( = (4 + 2a 2 ) /(x 2 + 4 + 2a 2 ). Since the lattice Green's function P&(C) is independent 
of the anisotropy parameter a, the critical point can be described by the parameter x 2 = 
2(2 + a 2 )(l — Cc)/Cc with ( c being independent of a, while T* oc 1/(2 + a 2 ). Then, utilizing 
Eq.fllo]) the critical density p c = x 2 T*/An is also independent of a. Similar arguments can 
now be applied for the analysis of the tricritical point. The thermodynamic behavior of 
the tetragonal body-centered lattice electrolytes is then different from that of the simple 
tetragonal lattice ionic systems. Since critical and tricritical temperature decay at the same 
rate as functions of the anisotropy parameter a, while p* and p* tri are constant, the possible 
low-density gas-liquid coexistence is suppressed by order-disorder phase transitions with a 
tricritical point at any degree of anisotropy, as shown in Fig.4. 



IV. DISCUSSION. 



Our analysis of simple tetragonal lattice models of electrolytes based on the DH ap- 
proach indicates that, similarly to isotropic sc ionic lattice systems, at weak anisotropies 
(0.385 < a < 2.113) the possible low-density gas-liquid phase separation is metastable, and 
the sublattice ordering is always thermodynamically more favorable. However, for strongly 
anisotropic simple lattices (a < 0.385 or a > 2.113) the gas-liquid coexistence reappears 
and the phase diagram becomes more complex, with critical, tricritical and triple points, as 
shown in Fig. 3. The explanation of this phenomenon is the following. For weakly anisotropic 
lattices the ordering of ions of opposite signs on different sublattices decreases significantly 
the total free energy, and order-disorder phase transitions with a tricritical point determines 
the phase diagram of the system. For strongly anisotropic lattices this ordering is less signifi- 
cant at low densities and the gas-liquid phase separation is restored. Another way of looking 
at this phenomenon, as was discussed above, is the analogy between lattice stretching and 
lowering of the space dimensionality. As was shown before,^ at low dimensions the critical 
temperature is increasing and thus gas-liquid coexistence occurs again at low densities. 

However, the thermodynamic behavior of Coulomb systems on body-centered tetragonal 
lattices is very different. At all degrees of anisotropy, the sublattice ordering is thermody- 
namically more stable, and gas-liquid phase separation is always suppressed. This is the 
result of the special symmetry of body-centered lattices. 

Our pure Debye-Hiickel treatment of anisotropic lattice electrolytes assumed that there 
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are only free ions and empty sites in the system. However, at low temperatures the formation 
of strongly bound neutral dimers, or Bjerrum pairs,lllll30 is a highly favorable process. Such 
pairing can be viewed as a reversible chemical reaction, and this process can be treated 
in a systematic way. Another important contribution to the free energy of electrolytes is 
the ion-dipole solvation energyB0'0 The formation of ion pairs and their interactions with 
single ions have a strong effect on the thermodynamic properties and phase diagrams of 
electrolytes.@0 Our theoretical method can be extended to include these effects, although, 
due to anisotropy of the lattice system, there will be more than one type of bound neutral 
dimers. Based on the comparison with the continuum and isotropic lattice electrolytesJliE 
we predict that, if we take these effects into account, both critical and tricritical temperatures 
will decrease, while critical and tricritical densities will increase. However, this means that 
our qualitative conclusions on thermodynamics and phase diagrams of anisotropic cubic 
lattices will not change. 

It is interesting to note that similar phase diagrams of lattice electrolytes have been 
obtained by Ciach and Stell.i They considered a mean-field theory of electrolytes with 
single-ion lattice potential on the isotropic sc lattice and additional short-range interactions. 
Note that this treatment neglects the cooperative screening, which is thermodynamically 
important and which is included in our DH-based theory. However, the origins of similar 
complex phase diagrams in both models are different. In the model of Ciach and Stell the 
gas-liquid phase coexistence is driven by short-range interactions, while in our model the 
lattice stretching makes the sublattice ordering less thermodynamically favorable at low 
densities and, as a result, gas-liquid phase separation is restored. 



V. CONCLUSIONS 



We have extended the Debye-Hiickel method to treat three-dimensional anisotropic lat- 
tice models of electrolytes. Phase diagrams for different degrees of anisotropy have been 
obtained. For weakly anisotropic simple tetragonal lattices, the order-disorder phase tran- 
sitions with a tricritical point suppress the possibility of low-density gas-liquid phase tran- 
sitions. However, for strongly anisotropic lattices gas-liquid phase coexistence is restored. 
Thus the lattice anisotropy for simple tetragonal lattices mimics the lowering of the spatial 
dimensionality. However, the thermodynamics of the body-centered tetragonal lattice ionic 
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systems is very different. There is no gas-liquid separation and the phase diagram has only 
order-disorder phase transitions for all degrees of anisotropy. This is the consequence of the 
special symmetry of body-centered lattices. 

The relevance of our results for understanding the thermodynamics of real ionic fluids 
remains unclear. However, our method may be more useful for description of the thermo- 
dynamics of real ionic crystals with defects.il Furthermore, numerical simulations of the 
anisotropic lattice ionic models with lattice Coulombic potentials are clearly needed in order 
to check the validity of our theoretical predictions. 
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Figure Captions 

Fig.l. Gas- liquid coexistence curves for simple tetragonal cubic lattices predicted by pure 
DH theory for different values of the lattice anisotropy parameter a. 

Fig. 2. Critical parameters as a function of degree of anisotropy: (a) critical and tricritical 
temperatures; (b) critical density; (c) tricritical density. 

Fig. 3. Phase diagrams of electrolytes on simple tetragonal lattices with sublattice or- 
dering for different degrees of anisotropy: (a) for a — 1 and a = 0.1; (b) for a — 15. 
Dashed lines show the metastable gas-liquid coexistence curves predicted by pure DH theory. 

Fig. 4. Phase diagrams for ionic systems on body-centered tetragonal lattices with sublattice 
ordering. The gas-liquid coexistence curves predicted by pure DH theory are shown by 
dashed lines. 
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Fig.l 
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Fig. 3a 
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Fig.3b 
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Fig.4 



